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ABSTRACT 

We study the interaction between gas and dust particles in a protoplanetary disk, comparing analytical and 
numerical results. We first calculate analytically the trajectories of individual particles undergoing gas drag in the 
disk, in the asymptotic cases of very small particles (Epstein regime) and very large particles (Stokes regime). 
Using a Boltzmann averaging method, we then infer their collective behavior. We compare the results of this 
analytical formulation against numerical computations of a large number of particles. Using successive moments 
of the Boltzmann equation, we derive the equivalent fluid equations for the average motion of the particles; these 
are intrinsically different in the Epstein and Stokes regimes. We are also able to study analytically the temporal 
evolution of a collection of particles with a given initial size-distribution provided collisions are ignored. 
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1. INTRODUCTION 

In an attempt to account for the coplanar nature of the orbits 
of all known solar-system planets, Laplace (1796) postulated 
that they were formed in a common disk around the protosun. 
Today, the detection of protostellar disks around most young 
T-Tauri stars (Prosser et al. 1994) is a strong evidence that the 
Laplace nebula hypothesis is universally applicable. The recent 
discovery of planets around at least 10% of nearby solar-type 
stars (Marcy et al. 2000) suggests that their formation may be 
a robust process. 

Conventional cosmogonical scenarios are based on the as- 
sumption that heavy elements in gas-phase condensed to form 
grains which then coagulated into planetesimals and grew into 
protoplanetary cores which can accrete, at least in some regions 
of the disk, massive gaseous envelopes around themselves (Pol- 
lack et al. 1996). The coexistence of gas and solid ice has been 
detected in some protostellar disks (Thi et al. 2002). In fact, 
protostellar disks are most conspicuous in their continuum ra- 
diation associated with the re-processing of stellar light by the 
grains (Adams, Lada, & Shu 1987). The apparent wavelength 
dependence in the thickness of the disk dust layer has been 
interpreted as evidence of grain growth (Throop et al. 2001, 
D'Alessio et al, 2001, Clarke et al, 2003) and settling (Shup- 
mgetal, 2003). 

The /im-to-cm continuum radiation signatures of the dust are 
observed to fade on the timescale of a few Myr (Beckwith, 
1999, Haisch et al 2001), signaling the depletion of grains in 
this size range. This suggests that heavy elements initially con- 
tained in this size range are either evaporated, ejected to large 
distance, accreted onto the host stars, or have coagulated into 
larger particles. The first possibility is constrained by the con- 
current decline in the CO-gas (Zuckerman et al. 1995) whereas 
the last possibility is directly relevant to the process of planet 
formation. 

Theoretical analysis suggests a very strong constraint on the 
growth of /im-size grains into km-size planetesimals. Indeed, 
the orbital evolution of the particles is determined by both the 
gravity of the central star and the drag of the disk gas. In the 
absence of turbulence, the disk gas attains a dynamical equilib- 
rium between gravity, pressure, and centrifugal forces with zero 
velocity in both radial and normal-to-the-disk directions and a 



slightly sub-Keplerian velocity in the azimuthal direction. Par- 
ticles in the disk undergo both sedimentation toward the mid- 
plane and inward drift in the radial direction (Whipple 1972, 
Weidenschilling 1977). In a minimum mass nebula (Hayashi et 
al. 1985), the resulting orbital decay timescale at 1 AU (for in- 
stance) is smallest for m-size particles (Adachi et al. 1976), and 
is then less than about 10 2 yr. Unless the growth of planetes- 
imals across this "most vulnerable size" can occur faster than 
their orbital decay, there would be no residual planetesimals left 
to provide the building blocks of planets. 

One possible channel of rapid grain growth is through sedi- 
mentation into a sufficiently thin, gravitationally unstable disk 
(Goldreich & Ward 1973). The critical thickness for gravita- 
tional instability of such disks is less than ~ 10~ 5 of their radii 
and the characteristic size of the resulting fragment is ~ a few 
km. However, even a modest amount of turbulence can pro- 
vide adequate stirring to prevent the sedimentation of grains 
into such a thin unstable layer (Weidenschilling 1984, Supulver 
& Lin 2000). Though turbulence is likely to occur in a magne- 
tized disk (Balbus & Hawley, 1990) through magneto-rotational 
instability, this mechanism could well fail in regions of the disk 
where the ionization fraction is too small. In these regions only, 
the following alternative mechanism for turbulence has been 
proposed. 

In a laminar disk, the sedimentation of dust toward the disk's 
mid-plane leads to a local concentration of massive particles; 
these particles entrain the gas to a near-Keplerian velocity 
through drag, thereby introducing a shear layer between the 
dust-dominated mid-plane and the rest of the disk gas (Wei- 
denschilling & Cuzzi 1993). Such a flow pattern in the disk has 
the potential to cause the onset of a shearing instability (Sekiya 
1998, Youdin & Shu 2002). However, the stability analysis 
used by these authors for such flow is based on a single-fluid 
approximation in which the dust particles are assumed to be 
well-coupled to the gas. Since the concentration of the dust 
particles not only causes the shear but also a stabilizing density 
stratification, the flow of dust and gas should be treated sepa- 
rately. In a companion paper (Garaud et al. in preparation), we 
will carry out a two-component stability analysis of the disk's 
dust layer. 

Such a study is greatly simplified by the treatment of the par- 
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tides as a separate fluid rather than a collection of particles. It 
is with this goal in mind that we now present a system of aver- 
aged equations for the evolution of a collection of dust particles 
in the form of moments of the Boltzmann equation. This pre- 
scription could also in principle be applied for the studies of 
dust particles' evolution due to coagulation, sublimation, con- 
densation (Supulver & Lin 2000) and under their interaction 
with embedded planets (Wyatt et al. 1999) and stellar radiation 
(Takeuchi & Artymowicz 2001, Klahr & Lin 2001, Takeuchi & 
Lin 2002). 

For the present calculation, we assume the particles are col- 
lisionless and indestructible spheres in a gaseous laminar disk 
with no embedded planets and negligible stellar radiation. In 
this paper, we also neglect the feedback of the particles' drag 
on the motion of the gas. In §2, we recall the general gas drag 
laws and their effects on particle trajectories in a protoplane- 
tary accretion disk. In §3, and §4 we solve for individual par- 
ticles' orbits, and derive the complementary set of dynamical 
equations for a collection of particles in the form of low-order 
moments of the Boltzmann equation (Boltzmann, 1872). In §3, 
we focus only on small particles, which have sizes smaller than 
the mean-free-path of the gas molecules, whereas in §4, we de- 
velop analogous equations for particles with sizes larger than 
the mean-free-path of gas. These initial calculations are purely 
analytic; in order to perform them, we need to assume separa- 
tion of the radial and normal-to-the-disk variables, as well as a 
constant background gas density. In §5, we check the validity 
of these assumptions with a complete numerical calculation of 
the particles' orbits in a simple-to-use model for the gas nebula. 
We compare both the results for individual orbits and for the 
moments (typically, the average velocity and average velocity 
dispersion) obtained from the Boltzmann averaging procedure. 
We also consider the spatial and temporal evolution of a given 
size-distribution of particles in §6 and show that particles with 
a narrow selected size-range sediment most rapidly toward the 
mid-plane. Finally in §7, we summarize our results and discuss 
their implications. 

2. GAS DRAG LAWS AND PARTICLE MOTION 

In this section, we briefly recall the effect of gas drag and 
gravitational forces on individual particle trajectories. 

2.1. Drag force description 

In dust-particle number-density regimes where particle colli- 
sions are negligible, the trajectory of a particle is given by 

r"=-V$-— F d , (1) 

m p 

where q' denote derivatives of q with respect to time, r is the po- 
sition vector of the particle, $ is the externally imposed gravita- 
tional potential, m p is the mass of the particle and Fd is the drag 
force between the particle and the gas. The amplitude of the 
drag force depends on the size of the particle s compared to the 
mean free path of the gas A and one typically distinguishes two 
regimes, s <C A (Epstein regime) and i>A (Stokes regime). 

In the small particle limit, i.e. the Epstein regime, the drag is 
caused by the thermal agitation of the gas and is proportional to 
the velocity of the particle relative to the gas: 

F d = m p — -(r-v g ) whens< A, (2) 

where v g is the gas velocity, p is the local density of the gas, c 
is the local sound speed and p s is the solid (i.e. internal) den- 
sity of the particle. In a standard solar nebula, this formula is 
typically valid up to decimeter-sized particles at 1 AU. 



In the other limit, large particles see the gas as a fluid, and 
experience a drag force through the laminar or turbulent wake 
that they create as they move through the gas. This is the Stokes 
regime. Whipple (1972) reported that the drag force on a sphere 
is 

F d = m p — ^— ^|Av|(r-v g ) wheni> A, (3) 

Ps s 

where | Av| is the norm of r-v g and Re is the particle Reynolds 
number of the flow (Re = 2s|Av|/V where v = Ac/3 is the 
molecular viscosity of the gas). Experimental results suggest 
that the drag coefficient C(Re) varies like 

C = 9Re~ l for Re < 1 , 

C = 9Re~ a - 6 for 1 < Re < 800 , 

C = 0.165 for 800 <Re. (4) 

This prescription ensures a smooth transition between the Ep- 
stein regime and the Stokes regime for intermediate particle 
sizes (s ~ A). Note that in the Stokes regime (for Re — > oo), the 
amplitude of the drag force depends on the total velocity of the 
particle with respect to the gas, and not simply its component 
parallel to the particle's velocity. This expresses the fact that 
a strong wake caused by the particles motion in one direction 
also affects even the slightest motion in another perpendicular 
direction. 

2.2. Separation of vertical and radial motions 

As particles undergo gas drag, they lose angular momen- 
tum which causes them to drift slowly inwards. For particles 
strongly coupled to the gas (i.e. very small particles), this drift 
is slowed down by the support provided by the gas pressure it- 
self. For weakly coupled particles (i.e. very large particles), the 
angular-momentum loss is very small and the drift is equally 
slow (see the next section). There exists an intermediate regime, 
however, where orbital decay can occur very rapidly. Weiden- 
schilling (1977) carried out a first quantitative study of the ef- 
fect of the gas drag on particles of various sizes, and deter- 
mined the inward drift velocity as a function of particle size. 
He found that the maximal drift occurs for particles for which 
the typical stopping timescale f, = |(r'-v g )|/|Fd| is of order of 
the orbital timescale , and decays very quickly for particles 
much larger, or much smaller than that. As a result, unless the 
particle is of intermediate size, there will be a marked sepa- 
ration between the dynamical timescale and the orbital decay 
timescale. 

Moreover, there is also a marked separation between radial 
and vertical length scales; indeed, numerical integration of the 
orbits of particles in a gaseous disk shows that the typical radial 
drift velocity is of the same order (in the case of small parti- 
cles) or much smaller (in the case of large particles) than the 
vertical settling velocity (see Section 5.2). Hence in a thin disk 
the particles settle to the mid-plane with very little radial excur- 
sion, then undergo radial decay within a very thin dust layer. 
For this reason it is justified to assume that the vertical motion 
of particles is more or less independent of their radial motion 
(although not of their radial position), and that their radial mo- 
tion occurs mostly at z = 0. Since Weidenschilling (1977) has 
already derived an analytical description of the particles motion 
in the radial direction, we shall concentrate on the problem of 
vertical settling. A short derivation of the case of radial and az- 
imuthal motion of the particles in a gas disk is presented in the 
Appendix for completeness. 
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3. THE COLLECTIVE EVOLUTION OF SMALL PARTICLES 

3.1. Vertical motion of small particles in the Epstein regime 

From the previous sections, we infer the vertical equation of 
motion of a small dust-particle in the Epstein regime: 



z"=- 



p c 



(5) 



dz p s s 

where z is the height above the disk. In a central potential, the 
gravitational force close to the mid-plane can be rewritten as 
— £l\(r)z where Vl^(r) is the Keplerian angular velocity at the 
radius considered. We re-normalize the time variable to the or- 
bital timescale t = Ojj at the radial position considered, and 
the distances to Astronomical Units (AU); we deduce the non- 
dimensional equation 

z = -z-pz, (6) 
where p = il^(r)(p(r)/ p s )(c(r)/s). In terms of this normaliza- 
tion, the magnitude of p in this regime is typically much larger 
than unity, which means that the stopping time 1 / p is much 
smaller than the orbital time. The solution to equation (6) is 
straightforward, and the trajectory of a particle initially posi- 
tioned at height z; with initial velocity w, is 



Z p (p.,t;Zi,Wi)= ^ 



Iwi + pzt 
z i +— F == | exp 




-p+^fp 1 



-//- 



Vp 2 



V Vp^J V 2 

Note that for p ^> 1 this expression simplifies to 



z ? (p,t;zi,Wi)= Zi+— e 



Wi 



-tit 



•(7) 



(8) 



P J P 

The second term of this expression describes the rapid deceler- 
ation of the particle by the gas drag, the first represents a slow 
settling toward the mid-plane. Let us recast this equation into 

z p (p,t;zi,Wi) = a(p,t)zi + P(p,t)wi , (9) 
which defines the functions a and (5 uniquely. The instanta- 
neous vertical velocity of these particles is then given by 

w p (fi,t;zi,Wi) = z p = dzi + Pwi . (10) 
Note that for p ~ 1 (which describes the intermediate param- 
eter range between the Epstein regime and the Stokes regime) 
the particle motion contains an oscillatory component decaying 
exponentially on timescale p/2. 

3.2. Boltzmann description and continuum equations 

Given a set of particles at initial positions z, and with ini- 
tial velocities w,, the Boltzmann distribution function of these 
particles is 

f(z,w,t) = ^2miS(z-z p (p,P,Zi,Wi))S(w-w p (p,t;zi,Wi)) . 

(11) 

One can also take the continuum limit of this description in the 
case of a very large number of particles: 

f(z,w,t) = J Azi J dwi pi(zi)g(wi,Zi)5(z-z p (p,t;zi,Wi)) 

■5(w-w ? (p,t;zi,Wi)) , (12) 
where p, is the initial particle density distribution, and g is the 
initial velocity distribution of particles. The freedom in the rela- 
tive choices of the initial distribution functions p, and g is lifted 
by requiring that for all z, 



/ 



dwi g(wi,Zi)= 1 . 



(13) 



The mass density of the particles is obtained by integrating / 
over all possible velocities: 

Pp(z,t) = J dwf(z,w,t) 

= J dzi J dwiPi(zi)g(Wi,Zi)5(z-z v (p,t;zi,Wi)) . (14) 

Substitution of the equations for the particle trajectories into 
(14) and integration over all possible initial positions yields 



P P (z,t) = - / dwipi [ — - — ) g [ wi, 



a 



(15) 



The first moment of the Boltzmann function is the average 
vertical velocity: 



Pp(z,t)w= J dwwf(z,w,t) (16) 
dz ; J dwiPi(zi)g(Wi ,Zi)WpS(z-Zp(p,t; z, , w,)) . 



The same manipulations yield 

_ fdwf ( z-0Wj \ ( z-/3wj \ ( . z-pwj A \ 
p p w = J — Pi{—)8{^—) [a—+Pw t ) 

(17) 

Let us assume for a short while that the distribution of the 
initial velocities is independent of height above the mid-plane. 
Then we can verify easily that the quantities determined above 
in equations (15) and (17) satisfy the standard continuity equa- 
tion 



(18) 



regardless of the form of the "trajectory functions" a and (3. 

However, the original collisionless Boltzmann equation 
(Boltzmann, 1872) describes the evolution of the distribution 
function / in phase space through 



dt dz dw 



(19) 



where T is a collision term that reproduces the interaction of the 
particles with themselves and with the surrounding medium. If 
we integrate the Boltzmann equation with respect to the veloc- 
ity space, we get 

^ + | ( p p W ) + /dww£ = /rdw. (20) 

Substitution of the equation of motion equation (6) into the last 
term of the left-hand-side provides an expression for the inter- 
action term T as a condition for the mass continuity equation to 
be satisfied: 



T = pf. 



(21) 



This result can be generalized in the case where g depends on 

Zi- 

Subtracting the mass continuity equation from the first mo- 
ment of the Boltzmann equation reveals the importance of the 
particle velocity correlations w 2 



dw dw 9$ 
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Pv^ + Pv w -fo=-Pv-Q^-PPv">-—SPv (J ^ W 



dz 
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where a 2 = w 2 — w 2 represents the particles' velocity dispersion 
function. This equation is very similar to that of a standard 
fluid, with the following caveats: (i) it contains an explicit drag 
term in the form of -pw, and (ii) the very last term, which usu- 
ally represents the pressure term in a standard fluid 1 , must be 
written explicitly as a function of known quantities. We expect 
this term to be null in the case of collisionless particles with 
negligible relative velocity with respect to the gas; in order to 
double check this conjecture, we now evaluate a 2 explicitly. 
Let us first evaluate w 2 : 



be used only in conjunction with a complete evaluation of the 
particle's orbit. However, one can assume that the largest typi- 
cal contrasts in velocity between particle and gas occurs in the 
vertical direction, where the particles oscillate across the mid- 
plane with an epicyclic frequency whereas the gas is more-or- 
less stationary (so that z' 2 3> (r6' - v g ) 2 and z a 3> r' 2 ). In that 
case, the expression for the vertical component of the drag force 
reduces to 



F d • e z = -m. 



' As 



z z 



(28) 



<,2- 



dwi ( z—(3wi 



z-PwA ( .z-fiwi ■ 
wt, a +/3w, 



a 



a 



(23) 

The behavior of a 2 at short times (i.e. for times shorter than 
the stopping time f s = l/p) is difficult to extract analytically 
while keeping the functions p-, and g unprescribed. However, 
for large p and large times one can simplify the expressions for 
the functions a and (3 to 

a = exp(-f / p) + 0(exp(-^f )) , 



(29) 



= - (exp(-*//i)-exp(-/i?)) 
A 



(24) 



In that case, the expressions for the first and second velocity 
moments simplify largely to 

z 

— + 0(sxp(-pt)) , 



Note that since the typical velocity of the particles across the 
mid-plane is of order of A17k where A is the typical height of 
particles above the disk, the asymptotic Stokes regime occurs 
for particles of size 

s _ Re c r v K 

a"Ta7 ' 

where Re c = 800. This condition corresponds for example to 10 
m size particles at r = 1 AU and A = 1 /10th of the disk scale 
height in a standard solar nebula. 

4.1. Particles' trajectories 

Using the same normalizations as in the Epstein regime, the 
particle's vertical equation of motion can then be written as 



w = - 



-~ z 



-z-p\z\z ■ 



(30) 



• + 0(exp(-pt)) , 



(25) 



where p = C(Re)(R / s)(p/ p s ) and R = 1 AU. In the asymptotic 
Stokes regime, C(Re) is simply a constant, and p is typically of 
order of unity, or smaller. 

This equation can be solved exactly through the introduction 
of the quantity w = z. On the positive branch (when w > 0) and 
negative branch (w < 0) the solutions are, respectively, 



wX=A exp(-2 P z) ~- + ^2 , 



w_ = B exp(2 pz) + - 



p 2p 2 
1 



w_ 



and one can show that provided g is an even function of velocity 

a 2 = <9(exp(-2/if)) . (26) 

This result is expected: any initial velocity dispersion is quickly 
damped out by the surrounding gas on the short stopping 
timescale f s . Note that in the case where p, and g are uniformly 
distributed, this expression is also valid for very short times. 
The long-term evolution of the collection of particles becomes 
a slow collective settling, on a timescale p. Therefore, in the 
Epstein case, one can approximate the particles as a fluid with 
a linear drag force and zero pressure. 

4. THE COLLECTIVE EVOLUTION OF LARGE PARTICLES 

The main drag effect of the gas on a large particle occurs 
through the formation of a wake behind the particle. If the ve- 
locity of the particle relative to the gas is small enough, the 
wake is laminar and the drag force is a linear function of the 
relative velocity. In that case, the particle trajectory is actually 
the same as that obtained in the case of the Epstein regime, and 
the results of the previous section apply. However, as the rel- 
ative velocity between the particle and the gas increases, the 
wake becomes turbulent and the amplitude of the drag force 
becomes a power law of the particle's Reynolds number. 

As a result, though the component of the drag force in any 
direction is still proportional to the relative velocities of the par- 
ticle and the gas in that same direction, its amplitude also de- 
pends on the total relative velocity of the particle with respect 
to the fluid. Hence the vertical component of the drag force is: 

F d • e z = -m p ^ ^ Jr* + z?Hr0'-v g W , (27) 
As s v 

if one assume the only important component of the gas veloc- 
ity is in the azimuthal direction. This formula can in principle 

1 this term is related to thermal elastic collisions of the gas particles in a standard fluid, and tends to depend only on the local density and temperature of the fluid 



(31) 



p 2p 2 

As the particle oscillates about the mid-plane, it follows one 
branch or the other. The constants A and B are determined by 
matching the first branch to the boundary conditions and the 
successive alternative ones to each other at the turning point 
(i.e. when w = 0). For instance, the trajectory of a particle start- 
ing from rest at z = Zo above the mid-plane satisfies the equation 

- Z j-^)^- 2 ^ )+L p + ^- (32) 

It accelerates toward the mid-plane, then decelerates on the 
other side until it reaches a turning point z\ which satisfies 

0= (---^- 1 ))exp(2pz 1 -2pz )+- + ^ 1 . (33) 
\ p 2p L ) p 2p L 

These solutions indicate the existence of two regimes. When 
the initial height zo is much larger than the scale height l/p 
the turning point zi is roughly equal to 1 /2p regardless of the 
height from which it originated. This pattern corresponds to a 
rapid stopping phase, when the particle quickly drops toward 
the mid-plane. When the amplitude of zo is reduced much be- 
low l/p, the turning point z\ is roughly equal to -zo, which 
corresponds to the oscillatory phase, when the particle follows 
an epicyclic motion about the mid-plane with a slowly decaying 
amplitude. 



5 



4.2. The very rapid stopping phase 

The stopping phase corresponds to the limit where z is much 
larger than the scale height 1 /p. Note that this situation may not 
always occur since in the Stokes regime, \i < 1 . This condition 
corresponds to particles which start at a considerable distance 
above or below the disk itself. Nonetheless, in this case the 
solutions simplify to: 

;z--^7tt, (34) 



2 Z 

w_ ~ — 



1/2 



(for the downward branch) and this differential equation can 
easily be solved with 

Zp(t)=(zl /2 -^) . (35) 



2/i»/2 



Note that the particle velocity in that case varies linearly with 
time: 



1/2 



t 



W ^ = - 2 {^) + 2~,> (36) 

so that the typical duration of the settling phase is f, = 2p}l 1 z^ 1 ■ 
Finally, we also note that for a particle starting from zq above 
the disk, for example, the lower turning point is more or less 
independent of zo and is roughly equal to z\ = -1 /2p. 

4.3. The slowly decaying oscillating phase 

The oscillating phase corresponds to the limit z«l//(. In 
that case, the equation for the trajectory of the particle (in the 
downward branch, for example) simplifies to 

w z _ = z 2 = 2z (zo - z) - ( 1 + 2 A1 z )(zo - z) 2 . (37) 
This first order differential equation can be solved with conven- 
tional substitutions to yield 



Zp(f) = Z0" 



1- 



1 + 2pzo 

The lower turning point is then simply 

2zo 



cos ((l+2 M z ) 1/2 f) 



Zl = zo - 



(38) 



(39) 



1 + 2^z 

Starting from that point, the upward branch is then obtained 
through a similar analysis and becomes 



z P (0 = zi - - 



Zl 



l-cos((l-2 A iz 1 ) 1/2 f 



(40) 



(41) 



1 - 2/xzi 

The following upper turning point is zi such that 

2zi 

Z2 = Zl ~ • 

1-2/izi 

These equations show that the amplitude 7(f) of the oscillation 
slowly decays with time, and can be approximated by 

3tt 
Apt + A 

where K is a constant which is determined in such a way as to 
satisfy the initial conditions. If the particle is released from rest 
from z = zo at t = (with z <C 1 / p), K = 3tt/zq. Moreover the 
frequency of the oscillation also varies with time, and converges 
to the epicyclic frequency. 

To recapitulate briefly, we find that after a brief stopping 
phase on a timescale proportional to /i 1//2 , the particle oscil- 
lates around the mid-plane with frequency which is close to the 
epicyclic frequency, and with an amplitude that decays slowly 
in time roughly as given by equation (42). These solutions are 
illustrated in Figure 1 . 



4.4. Continuum equations in the Stokes regime 

The true particle trajectory in the asymptotic Stokes regime 
cannot be expressed analytically in any simple way. As in the 
Epstein case, we will only look at long-term behavior. In that 
case, we can approximate the trajectory of a particle by 

z p (t,p;K i ,ip i ) = -f(p,t;K i )cos(t + ip i ) , (43) 

where Kj and iff are complicated functions of the initial position 
and velocity of the particle. For simplicity, we will consider the 
asymptotic limit for times 4/if >> Ki. In this limit, equation 
(43) with 7 = 3n/Apt is a good approximation, since the ampli- 
tude of the particle trajectory becomes independent of the point 
of release (see Fig. 1). Let us also assume for the purpose of 
this work that the phase simply has a given distribution function 
g(ip) so that Boltzmann function for the Stokes regime is 

f(z,w,t) = JdzoPiizo) J dipg(ip)S(z- z p (p,t;ip)) 

= ?, v [dipg(<p)5(z-z v (lJ,,t,<p)), (44) 

where E p is the column density of particles, and the normaliza- 
tion of g is chosen such that dtpg(tp) = 1 . The function g 
is periodic with period 2tt but the initial distribution of parti- 
cles can be described with < (p < 2ir. Using the well-known 
relation 

-1 



6(h(x)) = J25(x- Xi ) 



dh 

— (x = Xi) 

dx 



where h( Xi ) = , (45) 



we can determine the particle mass density profile, 

E / z 2 \" 1/2 

Ppi z,t)= -E W(7 -| z |)^i-_j 

■[s(cos _ V7)-0 + s(- coirl te/7)-0] ' < 46 > 
where Tt is a Heaviside function and the function cos" takes 
values in the interval [0,7r]. 

As in the Epstein case, we can compare the fluid equation for 
mass conservation to the first moment of the Boltzmann equa- 
tion and obtain an expression for the collision term in the Stokes 
regime: 

T = 2fi\w\f. (47) 

Similar operations with the second moment of the Boltzmann 
equation yields the equivalent of the fluid equation of motion, 

dw _dw 9$ /",,,, d 2x 
Pv ~df + PvW ~dz~ = ~ Pv ~dz~~J ^/kkdw-Oyr ) 
It is extremely tempting to write the drag term as 



(48) 



/ 



[if\w\w dw = np p \w\w . 



(49) 



However, a proper evaluation of the left-hand-side integral 
shows that this cannot be done formally, and the complete ex- 
pression for the drag force term should be kept. From equations 
(43) and (44), 



/ 



f\w\wdw=— ( 

7 V 7 2 



g(cos-\z/6)-t)[ £Z- 7 Jl-±- 



+ g(-cos l (z/6)-t) ^ + 7 
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FIG. 1 . — Left panel: Trajectories of two particles as integrated numerically from the differential equation (30). The particles were released from rest from heights 
lOr (solid line) and lOOr (dotted line), in the case where p, = 0.1. Note that the first lower turning point is always located around -l/2p and that the oscillation 
amplitude is independent of the point of release in the following oscillating period. However, depending on the initial condition, a phase difference can exist between 
two trajectories. Right panel: Trajectories of a particle released from rest from height 0. 1 above the disk, with /i = 0. 1 . The particle oscillates around the mid-plane 
with an amplitude that decreases algebraically with time; the dotted line was drawn according to the amplitude given by equation (42). 



In the case where g is uniformly distributed, this expression re- 
duces to 



so that in the asymptotic limit 



/ 



/|w|wdw 




- (51) 



whereas the simpler expression p v \w\w in equation (49), which 
is rewritten p P (z7/7)|z7/7| in that case, would be wrongly ap- 
plied. 

As before, we are interested in evaluating the velocity disper- 
sion a 2 . We need to calculate 

p27T 

^=r-= dipg(ip)6(z-^co&(if+t)){jcos(t + ip)-^&m(t + if)) . 
Jo 

(52) 

In the asymptotic limit 



PpW = W(t-JzJ) 
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-1/2 



^( v 1 -^) z-j)g(cos-\z/S)-t) 



-1/2 



z+j\g(-co S - l (z/8)-t) (53) 



= PLl z -H (l -\z\) 
\ 7 

• [g (cos" 1 (z/S)-t)-g(- cos" 1 (z/S)-t)] . 
Similarly, one can show that 

p p vvvv = S p / dip g(tp)6(z-~/cos(ip + t)) 
Jo 

■(jcos(t + (p)-^sm(t + ip)) 2 , (54) 



p v ww = p p 



2 -2 
z 7 2 
— 5" + (7 - 



z 2 ) 



-2£ p W(7-|z|)7Z 



• [g( C os-\z/S)-t)-g(-cos-\z/S)-t)\ . (55) 

In the case where the initial distribution function g is uni- 
form, we can simplify these expressions greatly and show that 
for \z\ < 7 (which corresponds to the thickness of the dust layer) 
then 

a 2 = ( 1 2 -z 2 )U{ 1 -\z\). (56) 

This expression means that within the dust disk, the particle ve- 
locity dispersion remains important at all times. The velocity 
dispersion mimics the effect of a pressure term which effec- 
tively slows down the settling. In this simplified case where 
g(ip) was taken to be a uniform distribution, it is actually possi- 
ble to relate a to intrinsic large-scale properties of the system. 
However in the more general case where g(ip) is not a uniform 
distribution there exists no simple relationship between a 2 and 
the local density of the gas as there would normally be in a stan- 
dard fluid. Instead, this mock-pressure terms would depend in 
a complicated manner on the initial configurations of the parti- 
cles in phase-space. This prevents any further progress, at this 
stage, in the use of a fluid description for these large particles. 

5. NUMERICAL CALCULATIONS 

With this simplified analytical approach, we have explored 
what effects the gas drag may have on particles of various sizes. 
Three principal approximations were performed: 

• we neglected the radial motion of the particles, and 
therefore implicitly assumed that it is possible to per- 
form a separation of the variables and of the equations 
of motion into individual components; 
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• we neglected the vertical variation of p, and used a sim- 
plified expression for the gravity, 

• we neglected the contribution of the azimuthal motion 
of the particles through the gas in the calculation of the 
gas drag in the Stokes regime. This approximation ef- 
fectively underestimates the gas drag. 

In order to check the validity of these approximations, we 
begin by integrating the orbits of particles in a realistic gaseous 
protoplanetary nebula, with the complete expression for the gas 
drag force. We compare them to the corresponding analytical 
expressions. We then release various single-size sets of 10,000 
particles in the disk from a small region and follow their evo- 
lution in time. At a given time t we then perform the Boltz- 
mann averaging procedure and compare the results with those 
obtained in Sections 3 and 4. 

5.1. Numerical Integration: method and model parameters 

The equation of motion of a particle in a gaseous accretion 
disk is given by equation (1). Expanding this equation into its 
components in a cylindrical coordinate system (r,9,z) yields 

r » = r9 >2 ^__!l.e r 

(r 2 + z 2 ) 3 / 2 m p 

e" = -2-e'-^--e e , 



GMz F d „ 
e z , 



(r 2 + z 2)3/2 mp ««• < 57 ) 

where for small particles, we use the Epstein drag force expres- 
sion 



^± = ^ C -(r'e r +(r9'-v g )ee + z%), 
m p p s s 



(58) 



and for very large particles, we use the Stokes drag force 

— = - Jr' 2 + (rB'- v g ) 2 + z' 2 (r'e r + if& - v g )e e + z'e z ) . 
m p p s s V 

(59) 

We normalize these expressions using R = 1AU = 1.5 x 
10 13 cm as unit distance, and Q^(R) = 2.0 x 10" 7 s as a unit time. 
This yields 

^ =re ' 2 - (r 2 + ^2 )3 /2 -^^)^ 

f a , . ( a v g (r,z)\ 



e '- 2 f-^r-^m) ■ 

f= ~ (r 2 + ^)3/2 ~^ r ' Z)Z "' (60) 

where v K (R) is the linear azimuthal Keplerian velocity at R = 
\AU and where in the Epstein regime, 

p(r,z) c(r,z) 



Pv{r,z) = 
and in the Stokes regime 



p s stt K (R) ' 



Vs(r,z)= C(Re)- 



(61) 



(62) 



The numerical integration the particles' trajectories in this 
model is performed using a fourth order Range-Kutta integra- 
tor. 

The quantities related to the disk structure (as the gas density 
p, the gas velocity v g and the sound speed c) are derived from 



the protoplanetary nebula model of Supulver & Lin (2000). In 
this model, the gas follows a polytropic equation of state with 



p=Kp 2 andc= ^f{Kp , (63) 

where p is the gas pressure, with K = 6.9 x 10 20 dyn cm 4 g" 2 , 
7 = 1 .4 for an ideal diatomic gas (composed of pure molecular 
hydrogen for example) and 



p(r,z) = 8.5 xl0" 12 (- J (l-z 2 /H 2 (r)) gem" 3 , 



(64) 



where r is in AU and H 2 (r) = 34 x l0- 12 K(r/6T 3 / 4 fl K (rT 2 R~ 2 
is the square of the disk scale height (in the dimensionless 
units). Finally, the dimensionless gas velocity is 



v g (r) v K (0 1 dp 



1 



v K (R) v K (R) 2p 8r v K (R)Q K (r) 



(65) 



_ r _i /2 8.5x10 gg 3 3/4^-1/4 _ r -l/2_ -1/4 

GM 4 ' 
if we ignore the variation of the gas density scale height with 
radius. This expression defines e uniquely; with the numerical 
values of K, G, given above and using M = 2 x 10 33 g we have 
e = 1.8 x 10~ 3 . This expression describes how gas pressure re- 
duces the effective gravity on the gas; this results in the slightly 
sub-Keplerian character of the gas velocity as found in equation 
(66). 

5.2. Numerical integration: particles trajectories 

5.2.1. Small particles 

Figure 2 shows the trajectories of small particles (0.1mm to 
lm) released at 1AU from height 0.01 AU, as calculated numer- 
ically using the procedure described in Section 5.1. It compares 
the results to those obtained analytically assuming separation of 
variables and constant p which are summarized in equation (7). 
The analytical fit is indistinguishable from the complete numer- 
ical solution, despite these approximations. 

5.2.2. Large particles 

Figure 3 shows the numerically integrated orbit of a 10m 
sized particle in a protoplanetary nebula, when released from 
rest at radius 1 AU and height 0.01 above the mid-plane. It com- 
pares the numerical results to the theoretical predictions given 
by equation (43), and here again, despite the approximations, 
the analytical fit is excellent. Note that integration for much 
longer times (about 1000 orbits) reveals a small but growing 
discrepancy in the amplitude of the oscillation. This is due to 
the mis-representation of the analytical expression for the am- 
plitude of the drag force compared to its true value (see equa- 
tions (27) and (28)). This discrepancy is discussed in more de- 
tail later. 

5.3. Numerical integration: Boltzmann averaging and 
collective behavior 

In this section, we now follow the evolution of a collection of 
10,000 uniformly sized particles, in the Epstein regime (of size 
1mm) and in the Stokes regime (with size 10m). We release all 
particles from a small interval in radius, for 0.9 < r < 1.1. The 
initial radial positions are uniformly distributed in that interval. 
Other initial conditions depend on the regime studied, in order 
to be best able to perform the adequate comparisons between 
theory and numerical experiment. After a time f, the particles 
positions and velocities are gathered, and binned into regular 
intervals in radius and in height. The total number of particles, 
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FIG. 2. — Top figure: trajectories of small particles in the (r,z) plane, released from radius 1 AU and heights z = ±0.01 AU successively above and below the 
mid-plane (to avoid crowding in the figure). The markers mimic a Pointcarre map, i.e. are positioned at the points where the particle crosses the 6 = plane. 
Particles sizes are respectively: ■ = 0.1mm, A = 1mm, □ = 1cm, A = 10cm, X = lm. Bottom figure: trajectories of the same particles in the (z,t) plane, with 
the same points-style coding for the particle sizes; again, the points represent positions of intersection with the 9 = plane. The analytical trajectories as given by 
equation (7), using simply the value of p at the position the particle, are shown as dotted lines: they are virtually indistinguishable from the numerically integrated 
trajectories. 
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FIG. 3. — Top figure: trajectory of a 10 m particle in the (r,z) plane, released from radius 1 AU and height z = 0.01 AU above the mid-plane (to avoid crowding in 
the figure). The markers mimic a Pointcarre map, i.e. are positioned at the points where the particle crosses the 9 = plane. Bottom figure: trajectory of the same 
particle in the (z,t) plane; again, the points represent positions of intersection with the 9 = plane. The analytical trajectory as given by equation (43), using simply 
an average value of p in that region (p$ = 0.08) is shown as dotted lines: they are virtually indistinguishable from the numerically integrated trajectories. Note 
that one must take into account the fact that the normalization of the time- variable varies with radius in the simple analytical model. The dashed line represents the 
envelope of the oscillation as given by equation (42). 
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the average velocity and the average second velocity moment 
are then calculated for each bin according to the formulae: 

p p (r,z) = m p N(r,z) , 
w(r,z) = ^2w p /N(r,z), 



ww(r,z) = ^2w 2 p /N(r,z) , 



(66) 



where the sums are carried out over all particles found within 
the individual bin centered on the position (r,z). 

When comparing the numerical data to the analytical solu- 
tions of the previous sections, we are careful to apply the correct 
normalization to the analytical formulae: indeed the normaliza- 
tion of the time variable depends on radial position in the an- 
alytical solutions, whereas in the numerical computations, it is 
normalized to the orbital frequency at 1 AU. 

5.3.1. Epstein regime 

Using the normalizations adequate to the numerical calcula- 
tions, the theoretical solutions for the average particle velocity 
and dispersion becomes 



_ zr 
w 



-3/2 



+ 0(exp(-fi(r)r 3/2 t) , 



a 2 oc <9(exp(-2^0)r 3/2 f) . (67) 

Figure 4 shows the very short time evolution of a collection of 
10,000 mm-size particles released with a uniform distribution 
in heights in the interval [-0.01 , 0.01] AU and a uniform distri- 
bution in velocities with -l.d-3 < w> < 1 A- 3 (in units of the 
Keplerian velocity at r = 1). For this figure, the binning in the 
averaging process is fairly coarse (10 bins in each direction) in 
order to generate adequate statistics and to visualize the results 
more clearly. The velocity dispersion is indeed found to decay 
exponentially in accordance with equation (67). 

The trajectory of these particles are then followed up to time 
t = 100 (i.e. about 16 orbits) and their average velocity is com- 
puted with a finer binning (50 bins in each direction). The re- 
sults are presented in Figure 5, and compared to the analytical 
solutions in equation (67). Once again in this regime, the ana- 
lytical solutions are found to match the numerical results very 
precisely. 

5.3.2. Asymptotic Stokes regime 

In the asymptotic Stokes regime, particles oscillate across the 
mid-plane with an algebraically decaying amplitude. Theoret- 
ical analysis of the particles' trajectories and velocities shows 
the velocity dispersion should remain important at all times. 

In order to simplify the comparison between the analytical 
solutions in Section 4 and the full numerical results, we begin 
by releasing 10,000 lOm-size particles in the disk with similar 
initial conditions as the ones that were used in the analytical 
work. Though these initial conditions may seem perhaps artifi- 
cial, they are the only ones that actually provided any analytical 
solution. We specify the initial values 

Zi = y cos (<a) . 



37T 3tt 
-4fJ-o-jpCOs(ipi)- — sm(ipi) . 



(68) 



the point of release. We choose these initial values in order 
to simulate a set of particles released with a large range of ve- 
locities and heights above the mid-plane. The gravitationally 
bound particles converge toward the mid-plane during the set- 
tling phase, and continue to oscillate coherently across the mid- 
plane with a random phase but roughly all the same slowly de- 
caying amplitude. The simulation is then started at t = when 
all the particles have an oscillation amplitude 3n/K, the cor- 
responding velocity as given by equation (68) with a random 
phase. In that case, using the same normalization we adopted 
in our the numerical simulations, we find the theoretical solu- 
tions for the velocity dispersion to be 

Figure 6 compares the results of the numerical calculations to 
the theoretical solutions given by equation (69). For relatively 
short time after the onset of the calculation, the numerically 
obtained dispersion is indeed extremely well approximated by 
the analytical formula. However, over much longer times, a 
systematic shift emerges, in which the numerical values are 
slightly lower than the theoretical solutions. This slight over- 
estimate of the velocity dispersion comes from the fact that the 
analytical solutions underestimate the total drag force by ne- 
glecting the contributions from the particles' motion in the ra- 
dial and azimuthal direction. In order to verify this conjecture, 
another numerical experiment is carried out in which the drag 
force is artificially set to 

F d = m v ^-C(Re)-\z\(re r + (re-v g )tg + z^) . (70) 

Ps s 

In that case, the contribution of the drag force in the vertical 
direction is exactly that used in analytical calculation. As ex- 
pected, the analytical expression is now able to reproduce the 
experimental results much more accurately. 

Having gained some experience with the Boltzmann aver- 
aging methods, and having moreover established the adequacy 
of the analytical approximations to the particles' trajectories in 
the Stokes regime, we now attempt to use the results of the 
numerical simulations to find the evolution of the particles' 
dispersion in the case where the initial conditions are chosen 
randomly (rather than in such a way as to allow easy analyt- 
ical calculations). In the follow up numerical computation, 
10,000 lOm-size particles are released with a uniform initial 
height and vertical velocity distribution chosen in the intervals 
[-0.01,0.01] AU and [-0.01,0.01]v K (fl). The azimuthal veloc- 
ity of the particles is set to be the Keplerian velocity at the point 
of release. With these initial conditions, the initial velocity dis- 
persion should be uniform within the dust layer. 

Figure 7 shows the evolution with time of this dispersion pro- 
file at radius r = 1 AU, when the particles are binned in 20 inter- 
vals in radius and heights. After a brief adjustment phase (not 
shown in the figure), the dispersion profile seems to be well ap- 
proximated by a Gaussian distribution. We attempt to match 
such a Gaussian profile to the numerical data points, and find 
that the best fit is given by 



<? 2 (z,t)-- 



4/zf + 37r/A ' 
A 2 (f) ( 2z 2 \ 
"AW ' 



■ exp 



(71) 



where ipt is uniformly distributed in the interval [0,27r]. The 
azimuthal velocity of the particles is the Keplerian velocity at 



where Ao = 0.01. This formula was obtained from the follow- 
ing heuristic arguments: according to the formula for the parti- 
cles' vertical motion (see equation 43), the typical vertical ve- 
locity of the particles near the mid-plane is of the order of A(f ), 
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FIG. 4. — Left figure: Evolution of the velocity dispersion profile for a collection of mm-size particles. The style of the points used corresponds to those in 
the right-side figure. Right figure Evolution of the velocity dispersion of the particles located around r = 1. The points correspond to the results of the numerical 
experiment, and the solid line is the analytical prediction of equation (67), using a value of [m which is an average value for this quantity around r = 1. 




FIG. 5. — Contour plot of the average particle velocity for a collection of 10,000 mm-size particles. The solid contours as well as the gray-scale correspond to the 
numerical results, and the dotted lines correspond to the analytical predictions of (67), using the local value of /i(r, z) for each bin. Particles were released with a 
uniform distribution of heights in the interval [-0.01,0.01] AU and of radii in the interval [0.9,1.1] AU. One can easily observed the settling and (negligible) radial 
drift in this diagram, as well as the excellent matching between theory and this numerical experiment. 
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FIG. 6. — Temporal evolution of the velocity dispersion of a collection of lOm-size particles, measured at r = 1 AU . The successive curves (from top to bottom) 
correspond to times t = 0, t = 2000, t = 4000, t = 6000 and t = 8000 in units of the angular period at this radius. The following style coding has been used: (■) for the 
realistic numerical experiment and (A) for the artificially reduced drag force (as given by equation (70)). The dotted lines correspond to the analytical predictions of 
equation (69). 



where A is the maximum height of the particles above the mid- 
plane. Moreover, A decays with time due to gas drag according 
to equation (42), so that one can readily deduce the evolution of 
A(f) for a set of particles released at t = from typical heights 
ranging between and A . The numerical factors have been 
fitted to the data. The analytical fit is also shown in Figure 7, 
and is found to reproduce the evolution of a satisfactorily. 

Finally, recalling that the analytical expression for the parti- 
cles' vertical motion slightly underestimates the true drag force, 
we compare the analytical fit given by equation (71) to an ar- 
tificial simulation in which the drag force is given by equation 
(70). In principle, the fit should be much better in this case. 
However, we find that the evolution of the dispersion in that 
case is better fitted by 

CT 2( z ,0 = ^A 2 (0ex P (-jy , (72) 

whilst keeping the expression for A(f) the same (see Figure 7). 
Though this difference is small, it still suggests that any attempt 
to find an analytical fit to the collective motion of Stokes parti- 
cles should be carefully calibrated before being used in further 
simulations. 

To conclude this Section, it was found that the analytical 
evaluation of one particle's orbit was extremely useful toward 
the determination of the collective motion of a large number of 
particles, and that the resulting averaged equations fit the nu- 
merical experiments extremely well. We have also been able to 
determine a simple heuristic way of closing the moments equa- 
tions at low order in the case of particles in the Stokes regime, 
which emulates the effect of the particles dispersion by a mock 
pressure term that can be related to large-scale properties (and 
initial conditions) of the fluid. 

6. EVOLUTION OF THE SIZE DISTRIBUTION OF PARTICLES 
IN A GASEOUS DISK 



Having discussed the time-evolution of a large number of 
particles of equal size, we now consider a population of 
different-size particles. Again, the Boltzmann approach proves 
useful and simple to use. Let us assume for simplicity that for 
a given size, and at a given radius, the initial particles to gas 
mass ratio is independent of z. Hence the initial distribution of 
particles of a given size s at time t = is 

Pi (s,z) = po(s)exp(-z 2 /2H 2 ), (73) 

where H is the typical gas disk height; for an isothermal disk, 
H = c/ril^- Let us also assume (again for simplicity) that the 
initial vertical velocity of particles is null. This idealized ad 
hoc setup could correspond to a condensation front where, at 
time t = 0, both hydrogen and heavy elemental gas diffuses into 
a cool region of the disk. The condensation and sublimation 
time scales are generally much shorter than the dynamical time 
scales (Supulver & Lin 2000) so that, at least the small grains 
may form with a similar density distribution as the gas and with- 
out any initial in-falling motion. The coexistence of large and 
small particles well above the mid plane requires runaway co- 
agulation of initially steady, microscopic dust grains into grains 
up to a size s. This rapid in situ growth is possible in regions 
where the particles have large surface density and much greater 
than unity area filling factor because the particles' collision and 
growth time scales are short compared with their dynamical 
time scale. The particles then start settling from their initial 
positions, and interact with the gas mostly with an Epstein drag 
law. 

The particles density distribution at time t is then simply: 

p p (s,Z,t) = -Ptiz/a) = po(i)exp(-z 2 /2A 2 )exp(f/,u(s)) , (74) 

a 

where a is defined by equation (9), and A = Ha = 
H exp(-f / fi(s)) is the evolving dust layer thickness. This equa- 
tion was derived from (17) using the ansatz (73) and g(Wi,zd = 
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FIG. 7. — Left panel: Temporal evolution of the dispersion profile for a set of 10,000 lOm-size particles released with a uniform spatial distribution and a 
uniform vertical velocity distribution. The various sets of points joined by a solid line correspond to times t = 2000, t = 4000, t = 6000, / = 8000, / = 10000 and 
t = 12000 (from the top-curve to the bottom curve). The initial velocity dispersion profile is not shown in order to preserve a readable scale on this plot. But it is 
approximately uniform in the dust layer as expected . The dotted lines are the corresponding analytical fits to the dispersion profile as given by equation (71). Right 
panel: Temporal evolution of the dispersion profile for a set of 10,000 lOm-size particles released with a uniform spatial distribution and a uniform vertical velocity 
distribution, with the artificial drag force amplitude given by equation (70). The curves are plotted at the same time-intervals as in the left-side panel. The analytical 
fit, as prescribed by equation (72), is shown in dotted lines. 




S(wi). We can deduce from this result the number distribution 
of particles of size s 

n(s,z,t) = n (^)exp(-z 2 /2A 2 )exp(f/^)) 

= n (s)exp(-z 2 exp(2ts/n)/2H 2 )exp(ts/K) , (75) 

where we have defined the constant n such that \i = n/s in order 
to write explicitly this equation in terms of its dependence on 
the particle size. The initial particle size distribution function 
is taken to be that given by Hellyer (1970) and Mathis et al. 
(1977) 

n (s) oc s~ 3 5 . (76) 

At a given height above the disk, this distribution function 
peaks for small particle sizes (which corresponds to the initial 
peak in no(s)), but has another maximum, which evolves with 
time, at approximately 



regardless of the initial size distribution function of the parti- 
cles. This result is illustrated in Fig. 9. 

These results shown in Figures 8 and 9 suggest the following 
comments. Firstly, at a given height, successive fronts of parti- 
cles with decreasing sizes pass through the gas as they settle, to 
leave the region completely depleted of the intermediate sized 
particles. Secondly, at a given time, there is a definite strati- 
fication of the particles in terms of their respective sizes, with 
the larger particles concentrating tightly around the mid-plane. 
Note that this method is only valid for particles with sizes up to 
the order of the mean-free-path of the gas. For the very large 
particles, the Stokes drag law must be applied. In that case, 
we have shown in Section 4 that the large particles only settle 
algebraically toward the mid-plane. As a result, although their 



number density is small, the large particles remain important at 
all times and at all heights in the disk. 

This work shows that after a few settling times (which are of 
the order of a few thousands of orbits for cm-sized particles), 
the bulk of the disk will be depleted of intermediate sized par- 
ticles, which have all settled down to the mid-plane. There re- 
mains, at all heights, very small particles that have not had time 
to settle, and very large ones which are in a phase of oscillation 
around the mid-plane. 

7. SUMMARY AND DISCUSSIONS 

In this paper, we provide a set of analytic solutions to de- 
scribe the sedimentation of dust particles due to gas drag in a 
protostellar accretion disk. 

The main effect of the gas drag is to induce a loss of an- 
gular momentum and linear momentum for the particle, which 
in turn leads to orbital decay and vertical settling. Numerical 
integration of the particles' trajectories in a protoplanetary ac- 
cretion disk is straightforward if we assume that the gas flow 
is laminar and occurs mostly in the azimuthal direction. How- 
ever, analytical integration is also possible given a set of as- 
sumptions (which have been checked here to be reasonably well 
justified), including the separability of the radial and vertical 
motion, and locally uniform density. Accordingly, we have ex- 
tracted analytical formulae for the vertical motion of particles 
in a disk, both in the Epstein regime (for particles with size 
much smaller than the mean-free-path of the gas molecules) 
and in the Stokes regime (for particles much larger than the 
mean-free-path). Comparison of these analytical expressions 
with the complete numerically integrated orbits shows an ex- 
cellent match in both cases. 
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FIG. 8. — Evolution with time of the particle size distribution n(s) at fixed height z = 0.001 (i.e. 1% of the disk scaleheight). Note how the maximum of the 
distribution function satisfies the relation (77). 
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By using the analytical expression for the particles' orbits, 
we were able to extract their average behavior using a stan- 
dard Boltzmann analysis (i.e. local averaging, or rather, coarse- 
graining, of the particles' positions and velocities in phase 
space). The successive moments of the Boltzmann distribution 
function provide the evolution of the coarse-grained particle 
density, velocity and velocity dispersion. The successive mo- 
ments of the Boltzmann equation yield the equivalent of "fluid" 
equations for the collection of particles, and in particular a mass 
conservation equation and a momentum equation. 

Numerical experiments were performed in which 10,000 
single-size particles were released in the disk from a localized 
region of phase-space (typically uniformly distributed in space 
and with a uniform vertical velocity distribution). Using a sim- 
ilar coarse-graining method as the one used in the analytical 
analysis, we were able to extract the particle density profile, 
their average velocity and velocity dispersion from the numer- 
ical integration. We compared them to the analytical formulae 
obtained previously. The main results are summarized here. 

In the Epstein regime, there is an excellent agreement be- 
tween the analytic formulae and the results of the numerical 
calculation. Individual particles have their velocity with re- 
spect to the gas quickly damped out on a very short stopping 
timescale. Correspondingly, their collective behavior is coher- 
ent: the typical velocity dispersion is found to decay on the par- 
ticle stopping timescale, which is usually much smaller than the 
orbital timescale. For times longer than the stopping timescale, 
individual particles settle exponentially toward the mid-plane 
and the corresponding average momentum equation reduces to 
a simple standard fluid equation with a linear drag term, and 
null pressure. 

In the Stokes regime, there is a good agreement between the 
theory and the results of the numerical experiments. Individual 
particles are found to oscillate across the mid-plane with an al- 
gebraically decaying amplitude. Although this behavior is very 
well reproduced by the analytical work, there exists a small dis- 
crepancy between analytical formulae and the results of numer- 
ical computation in the form of a slight systematic overestimate 
of the particles' oscillation amplitude. This problem is due to 
an (unfortunately unavoidable) over-simplification of the am- 
plitude of the drag force in the analytical work, and results in 
similar discrepancies in the comparison of averaged quantities. 
Nonetheless, apart from this small identifiable and controllable 
error, the predictions of the analytical work match very well 
the numerical results. As expected, the particles' velocity dis- 
persion is found to remain large at all times within the dust 
layer, with a maximum near the mid-plane. This large disper- 
sion is directly related to the continuous oscillatory motion of 
the particles across the mid-plane. We found it extremely diffi- 
cult to obtain analytical predictions of the temporal evolution of 
the dispersion for anything but the simplest initial conditions - 
which were not necessarily always realistic. However, we were 
able to deduce from the numerical experiments a heuristic fit to 
the observed evolution of the particles dispersion, for more gen- 
eral initial conditions. This result has interesting consequences 
with respect to the possibility of closing the moments equations 
at low order, and describing accordingly the evolution of a col- 
lection of large particles with averaged fluid equations. We will 
test this theory in future work. 

Finally, by combining the average evolution equations for 
single-size particles obtained in the Epstein regime to a plau- 
sible size-distribution function for dust-grains in the ISM 



(Hellyer, 1970, Mathis et al., 1977), we showed that it is possi- 
ble to follow analytically the spatial and temporal evolution of 
such a size-distribution as the particles settle toward the mid- 
plane. This analytical calculation was done for collisionless 
particles, and the extension of this work to include collisions, 
coagulation and sublimation would probably require numerical 
analysis. Nonetheless, within this approximation we found that 
the sedimentation process results in a very strong segregation 
of particles according to their sizes, within a timescale of a few 
hundreds of years only. All intermediate-size particles quickly 
converge to the mid-plane, leaving behind only the very small 
particles, indirectly supported against the settling by gas pres- 
sure, and the very large particles, which keep oscillating across 
the mid-plane with a slowly decaying amplitude. This segre- 
gation process may well suggest that within the very thin dust 
layer at the mid-plane, one could approximate the population of 
dust particles by a single-size population. 

In future work, we shall apply the results obtained here to 
two essential purposes: 

• the description of the particles as a continuum fluid will 
be used to study the stability of the dust layer against 
shearing (Garaud et al, in preparation), in order to test 
whether gravitational instability can indeed take place in 
the dust layer. More generally, we intend to implement 
this description into existing 3D HydroDynamics codes 
for the computation of accretion disks, in order to obtain 
a fully self-consistent description of dust evolution and 
growth in these disks; 

• the temporal and spation evolution of the size- 
distribution of particles can be compared directly with 
observations of dust. We plan to construct models with 
which we can infer the extent of particle coagulation 
and infer the surface density distribution of the gas from 
high resolution, multi-wavelength (mostly in the sub- 
mm and mm range) maps of protostellar disks. 
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APPENDIX: DERIVATION OF THE PROPERTIES OF RADIAL 
MOTION OF PARTICLES IN A GASEOUS DISK 

Neglecting the vertical motion of the particles, the equations 
of motion in the radial and azimuthal directions in the Epstein 
regime are 

•o GM c p 

r-r9 2 = , sLf, (78) 

r s p s 

r6 + 2f0 = --—( re-vJr)) , (79) 

where v g (r) is the slightly sub-Keplerian azimuthal velocity of 
the gas 

fGM 1 1 dp n(r) 
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which defines ij(r) = -(1 /2p)dp/dr. As in Section 2.3, we use 
the normalizations [r] = R = 1 and [t] = fl^(R) to obtain 



f— r9 2 = — \ — fir , 



r9 + 2r9 = - 



r0- 



1 r](r) 



(81) 



(82) 



rV2 

It is not possible to find simple analytical solutions of the 
full system. However, we can assume that deviations from Ke- 
plerian motions are small, and that the particle velocity can be 
written as 

r = l + e and = 1 + 0. (83) 

The linearized system yields (assuming that the background 
quantities vary very little with r) 

e+fj,e-3e = 24>, (84) 



• ) + fJi( j) = -2e--fj,e-fj,-^- 



(85) 



Substituting the first into the second yields 
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_e+ M e+-(l + /x 2 )e = -/x^ r 



hence the local solution is 



i = e A "(acos2f + /7sin2f)-2 



i+/j 2 n| ' 



(86) 



(87) 



This solution represents the local difference between the parti- 
cle velocity and the Keplerian velocity, in units of the Keplerian 
velocity. As in the solution found in Section 4, the first term 
represents a very quick stopping of the particle on timescale 
1 /n and the second term represents a constant inward drift. The 
variation of this drift velocity with particle size is given by the 
function fi/(l + /j 2 ) which has a maximum for \x = 1 . This result 
reproduces the results found by Weidenschilling (1977). 
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